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The burst-suppression pattern is well recognized as a distinct feature of the mammalian 
electroencephalogram (EEG) waveform. Consisting of alternating periods of high ampli- 
tude oscillatory and isoelectric activity, it can be induced in health by deep anesthesia as 
well as being evoked by a range of pathophysiological processes that include coma and 
anoxia. While the electroencephalographic phenomenon and clinical implications of burst 
suppression have been studied extensively, the physiological mechanisms underlying its 
emergence remain unresolved and obscure. Because electroencephalographic bursting 
phenomenologically resembles the bursting observed in single neurons, it would be rea- 
sonable to assume that the theoretical insights developed to understand bursting at the 
cellular ("microscopic") level would enable insights into the dynamical genesis of bursting 
at the level of the whole brain ("macroscopic"). In general action potential bursting is the 
result of the interplay of two time scales: a fast time scale responsible for spiking, and 
a slow time scale that modulates such activity. We therefore hypothesize that such fast- 
slow systems dynamically underpin electroencephalographic bursting. Here we show that 
a well-known mean field dynamical model of the electroencephalogram, the Liley model, 
while unable to produce burst suppression unmodified, is able to give rise to a wide vari- 
ety of burst-like activity by the addition of one or more slow systems modulating model 
parameters speculated to be major "targets" for anesthetic action. The development of a 
physiologically plausible theoretical framework to account for burst suppression will lead to 
a more complete physiological understanding of the EEG and the mechanisms that serve 
to modify ongoing brain activity necessary for purposeful behavior and consciousness. 
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1. INTRODUCTION 

Prior to the development of the modern intensive care unit in 
the early 1960s, that featured intubation, artificial respiration, and 
comprehensive physiological monitoring, reports of the electroen- 
cephalographic pattern of burst suppression (BS) were confined 
to animal studies involving deep anesthesia and the occasional 
case of psychosurgery (Niedermeyer, 2009). Since then the burst- 
suppression pattern has become well recognized as a major diag- 
nostic feature of the EEG waveform that is encountered in a range 
of encephalopathic conditions, in addition to its appearance in 
health during deep anesthesia. Typically the BS pattern consists of 
bursts of high amplitude slow, sharp, or spiking electroencephalo- 
graphic activity separated by periods of electroencephalographic 
suppression (isoelectricity). The oscillatory features of the bursts, 
together with their duration and the duration of suppressed peri- 
ods show a high degree of variability (see Figure 1 for examples) 
that presumably reflects its myriad of initiating causes. First iden- 
tified during deep anesthesia with tribromoethanol in cats (Der- 
byshire et al., 1936), labeled burst-suppression pattern by Swank 
and Watson (1949) during barbiturate and ether anesthesia in 
dogs, it is now associated with cortical deafferentation (Henry 
and Scoville, 1952), cerebral anoxia and hypoxia, various types 
of intracortical lesions (Fischer- Williams and Cooper, 1963), deep 



coma, various infantile encephalopathies, the final stages of deteri- 
orated status epilepticus (Treiman et al., 1990), hypothermia, and 
high levels of many sedative and anesthetic agents (Schwartz et al, 
1989; Akrawi et al, 1996). 

Burst suppression in the absence of anesthesia is in general 
associated with a very poor prognosis. For example in neonates 
(Grigg-Damberger et al., 1989) the appearance of BS, even if tran- 
sient, is a portent of death or severe neurodevelopmental disability. 
In contrast, in adult populations while an anoxic/hypoxic BS pat- 
tern signals a serious pathophysiological event the outcome is not 
necessarily fatal and recovery with or without severe neurological 
damage is possible (Niedermeyer, 2009). Consistent with this are 
results of experimental work with EEG monitoring in rats reveal- 
ing that animals with greater rates of high amplitude bursts have a 
better survival and neurological outcome compared to those with 
lower rates of low amplitude bursts (Geocadin et al., 2002). 

While the electroencephalographic phenomenon and clinical 
implications of BS have been studied extensively (Brenner, 1985; 
Niedermeyer, 2009) the physiological mechanisms underlying its 
emergence remain in general unresolved and obscure. 

Burst suppression is typically thought to be spatially homoge- 
neous with burst onset and termination reported to occur near 
simultaneously across the entire scalp (Brenner, 1985; An et al., 
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FIGURE 1 | Example traces of electroencephalogram and 
electrocorticogram illustrating the heterogeneity of BS patterns. 

(A) Changes in neocortical electroencephalogram in the rat, recorded 
using dural surface electrodes, in response to a 5mg/kg/min thiopental 
infusion [figure reproduced with permission from Lukatch and Maclver 
(1996)]. (B) Electroencephalogram recorded in acute anoxia showing a 
clear burst-suppression pattern with grouped spikes [figure reproduced 
with permission from Hockaday et al. (1965)]. (C) BS pattern during 
closed loop target controlled propofol infusion at a target level of 



approximately 15/jg/ml (data courtesy of Professor Michel Struys, 
Groningen). Note the bursts consist of fast activity (>10 Hz) on a slow 
wave background. (D) Electrocorticogram obtained from an adult 
merino sheep during deep enflurane anesthesia, demonstrating high 
amplitude spikes interspersed with isoelectric periods of variable 
length [figure reproduced with permission fromVoss etal. (2006)]. (E) 
Electroencephalogram recorded from a 3-month-old infant suffering 
from infantile myoclonic encephalopathy [reproduced with permission 
from Niedermeyer (2005)]. 
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1996; Ching et al., 2012), indicating that low level subcortical 
mechanisms may be playing a decisive role. However arguing 
against this is the fact that this pattern persists subsequent to corti- 
cal deafferentation (Lukatch and Maclver, 1996), indicating that it 
probably represents an intrinsic, though physiologically abnormal, 
dynamical mode of cortex. Indeed the phenomenal resemblance of 
the patterns of BS to disorders of neuronal hyperexcitability sug- 
gests the involvement of similar physiological mechanisms. For 
example the bursting during burst suppression is often associated 
with myoclonic jerks resembling those seen during epileptic fits. 
Like generalized epileptiform activity, bursts are recorded simul- 
taneously at multiple electrode derivations, implying the wide 
synchronization of neuronal activity. 

At the cellular level a commonly reported finding is that hyper- 
polarization of the membrane potential of cortical neurons reli- 
ably precedes any overt electroencephalographic activity of BS 
(Steriade et al., 1994). Such hyperpolarization, which has been 
attributed to an increase in neuronal membrane potassium con- 
ductance (Steriade et al, 1994), has been hypothesized to play a 
major role in the induction of BS. This implied importance of 
inhibition in the genesis of BS is further supported by results 
involving rat neocortical brain slice micro-EEG preparations in 
which the application of a direct acting GABAyi agonist, musci- 
mol, readily induces BS. However contradicting this result is work 
reporting that inhibition is diminished during isoflurane-induced 
BS, in an in vivo feline preparation, as evidenced by increases in 
cortical neuronal input resistance and extracellular chloride con- 
centration (Ferron et al, 2009). Of course it may be that slow 
periodic modulations in inhibition, rather than singular increases 
or decreases in inhibition, underpin BS. In support of this view 
is the recent model of Ching et al. (2012), in which alterations 
in brain metabolism, due to the effects of hypoxia or anesthesia, 
parametrically regulate an activity dependent slow modulation of 
an adenosine triphosphate-gated potassium channel conductance 
to give rise to BS. However modulations in inhibitory activity 
alone may not be sufficient to account for BS and more consider- 
ation might need to be given to other mechanisms. For example 
Rroeger and Amzica (2007) present empirical evidence suggest- 
ing that modulations in excitatory synaptic efficiency, due to the 
progressive depletion of interstitial calcium during the periods of 
high amplitude electroencephalographic activity and its recovery 
during isoelectric periods, might account for BS. Consistent with 
this are reports involving laboratory slice preparations in which 
burst suppression induced by thiopental, propofol, and isoflurane 
is abolished by the application of glutamate receptor antagonists 
(Lukatch and Maclver, 1996). Whatever the pathophysiology of BS 
is it is reasonably clear that it is unlikely to be accounted for by a 
unitary physiological perturbation. That the physiological factors 
identified to date in BS all lead to a single well defined state sug- 
gests the possibility of an unifying dynamical mechanism. Thus 
the best hope for progress in understanding the phenomena of BS 
may be theoretical. 

How might we theoretically approach BS? The well studied 
dynamical mechanisms of action potential bursting (Izhikevich, 
2007) maybe able to provide vital insights into the mechanisms of 
bursting in the EEC In general the dynamical mechanisms under- 
lying bursting can be divided into two broad classes (i) fast-slow 



bursters in which there is a clear separation of the underlying 
time scales, with a fast system responsible for the fast spiking, and 
a slow system its slow modulation, and (ii) "hedgehog" bursters 
(Izhikevich, 2000) in which there is no clear separation of time 
scales. In terms of developing a theory of BS the former might 
represent the preferable starting point as the little empirical evi- 
dence that is available (Ching et al., 2012), at least in humans, 
suggests that alphoid activity, indicative of normal resting EEC, is 
preserved during the bursts of BS. Thus a theoretical starting point 
to understanding BS might be to consider the slow modulation of 
a dynamical system developed to describe the resting EEC 

One such dynamical system is the mesoscopic electrocortical 
model of Liley etal. (Liley etal, 1999,2002, 2011; Bojak and Liley, 
2005; Frascoli et al., 2011). This model is capable of accounting 
for a range of resting electroencephalographic phenomena that 
includes the alpha rhythm (Liley et al, 2002), the modulation of 
resting activity by sedative and anesthetic action (Bojak and Liley, 
2005) as well as the proconvulsant properties of the latter (Liley 
and Bojak, 2005), all within a physiologically plausible/admissible 
parameter space. This model is therefore well suited as a founda- 
tion from which to explore the physiological and dynamical genesis 
of BS. However, because in this model rhythmogenesis emerges 
from a strong coupling between cortical excitatory and inhibitory 
population activity, in its present form it has a restricted ability to 
exhibit BS through the parametric separation of time scales, either 
through the simulated actions of anesthetics or through other 
parametric routes. Here we show that BS can however emerge in 
this model by the addition of a slow system driven by one or more 
of the originally defined mean fields. We speculate that such a 
slow system represents a mathematical ansatz for the slow neuro- 
modulation of activity by a variety of intracortical, inter-cortical, 
and subcortical systems that include thalamus and the ascending 
neurotransmitter modulatory systems. 

2. MATERIALS AND METHODS 

2.1. MESOSCOPIC MEAN FIELD MODELING OF ELECTROCORTICAL 
ACTIVITY 

The electroencephalogram and electrocorticogram arise out of 
the cooperative activity of many thousands of neurons. A sin- 
gle electroencephalographic electrode records the synaptically 
induced currents of well over a 100,000 neurons (Nunez and Srini- 
vasan, 2005) and thus detailing each neurons contribution to this 
summed activity would appear superfluous. For this reason it is 
preferable instead to model the activity of populations of neurons. 
One general way of achieving this, in which known stochastic fluc- 
tuations can be included, is to dynamically evolve the probability 
distributions associated with the states of the neuronal ensemble. 
While in principle providing a rigorous way forward the formula- 
tion of such stochastic equations of motion entails a great deal of 
physiological uncertainty. For this and other reasons (Deco et al, 
2008) a more resolute path is to dynamically evolve some average 
quantity such as the mean soma membrane potential or the mean 
firing rate of some suitably defined neuronal ensemble. In this 
manner a mesoscopic level model can be developed which acts as 
a bridge between cellular (or microscopic) level activity and whole 
brain (or macroscopic) level behavior. While the current math- 
ematical approach for formulating the equations of motion for 
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FIGURE 2 | Schematic overview of the essential intracortical and 
cortico-cortical interactions between excitatory and inhibitory 
neuronal populations in the model of Liley et al. (Liley et al., 1999, 
2002, 2011; Bojak and Liley, 2005; Frascoli et al., 2011). Following 
conductance based approaches typically used to model networks of 
synaptically interacting networks of individual model neurons, excitatory 
and inhibitory neuronal populations are modeled as single passive 
resistive-capacitive compartments into which all synaptically induced 
postsynaptic currents flow. Functionally these populations are equivalent to 
the excitatory and inhibitory KO sets of Freemans K-set hierarchy (Freeman, 
1975). Cortical activity is then described by the mean soma membrane 
potentials of the spatially distributed excitatory (hj and inhibitory (/],) 
neuronal populations. The connection with physiological measurement is 
obtained through h„ which is assumed to be linearly related to the surface 
recorded electroencephalogram (Freeman, 1975; Nunez and Srinivasan, 
2005). Figure reproduced with permission from Frascoli et al. (2011). 



the activity of neuronal populations or masses, stems principally 
from the works of Wilson and Cowan (1972, 1973), Nunez (1974), 
Freeman (1975), and Amari (1975), they are not particularly suc- 
cessful in articulating the genesis of rhythmic activity in the EEC 
and its modulation by pharmacological agents, due to a range of 
mathematical simplifications that cannot be justified by an appeal 
to physiology. For this reason a range of biologically more refined 
neuronal populations models have been developed (Deco et al, 
2008) that form a more appropriate basis from which to model 
the physiological genesis of the electroencephalogram. The model 
of Liley et al. (Liley et al, 1999, 2002, 201 1; Bojak and Liley, 2005; 
Frascoli et al, 201 1 ) is one such model and will be the focus of our 
attempts to account for the phenomenon of BS. 

The model of Liley et al. is able to account for the major features 
of the mammalian electroencephalogram in the context of a para- 
meterization that accords with physiological measurement and 
experiment. Because it aims to provide a dynamical description 
of the electroencephalogram the model is spatially coarse-grained 
over the approximate extent of a cortical macrocolumn, reflecting 
the general physiological wisdom that such columniation signi- 
fies populations of neurons having similar functional behavior 
and anatomical connectivity. The multiple interactions between 
individual neurons are replaced by effective feed-forward and feed- 
back interactions between the mean activity (or mean fields) of 
excitatory and inhibitory populations of neurons over the spa- 
tial extent of the column. In this way mammalian neocortex is 
conceived as consisting of localized populations of excitatory and 
inhibitory neurons interacting by all possible combinations of 
feed-forward and feed-back connectivity interacting globally via 
long-range excitatory connections (Figure 2). 

Thus the response of the mean soma membrane potential 
h k (k = e, i) at position r on a two-dimensional cortical sheet, 
in response to induced post synaptic activity l\ k (/ = source, 
k = target population) is given by 

*k =h r k -h k (r,t)+J2 JTeq 771 fcfo t) (D 

l=e,i I 'k n k\ 

The postsynaptic response to a single pre-synaptic action 
potential (at f = 0) is modeled by the well-known synaptic alpha 
function of cable theory as r^y^t exp(l — yi k t)®(t) where 
is peak amplitude (occurring at t = f pe ak = 1/y/jfc) of the respec- 
tive excitatory (/= e) or inhibitory (/= i) postsynaptic potential 
(PSP), and 6(f) is the Heaviside step function. Thus we assume 
that the time course of the synaptically induced excitatory and 
inhibitory currents is described by a critically damped oscillator 
driven respectively by the mean rate of incoming excitatory and 
inhibitory axonal pulses: 

(J^ + Ylk^J Ilk(r, t) = exp(l)T, k y lk A lk (r, t), (2) 

with 

A ek (r, t) = N P ek S e [h e {r, f)] + 4>ek + Pek(r, 0. (3) 

and 

A tk (r, t) = Nf k Si[h(r, t)] + p ik (r, t), (4) 



where Ai k comprises the different sources of pre-synaptic spikes: 
NtlSi (input from local cortical neuronal populations), (p ek (input 
from long-range excitatory cortico-cortical fibers), and pi k (extra- 
cortical sources). While the present consensus is that extra-cortical 
sources (thalamo-cortical afferents) are purely excitatory in nature 
and thus p{ k = 0, we choose to retain these terms as when time 
independent they can be utilized to include the effects of tonic 
inhibition that are known to be induced by anesthetic action. The 
time courses of the synaptically induced currents, \ly\ k are taken 
to describe the time course of "fast" excitatory [l = e: a-amino-3- 
hydroxl-5-methyl-4-isoxazole-propionate (AMPA) and kainate] 
and inhibitory [l = i: y-amino-butyric-acid type A (GABAa)] 
neurotransmitter kinetics. Thus each type of PSP (excitatory, 
inhibitory) is described by two parameters Fi k , y\ k . However, as 
we will describe later, a parametrically more flexible description 
of the PSP is required to meaningfully model the effects of anes- 
thetics in which we can independently vary peak amplitude, rise 
(fpeak) and decay times. Mean neuronal population firing rates, 
Si, are assumed to instantaneous sigmoid functions of the mean 
soma membrane potential i.e., 

S,[h,(r, f)] = S ; ma 7{l + exp[V2(/i ; (r, t) - A/M)]} (5) 

The axonal pulses <p ek propagated by the exclusively excitatory 
long-range cortico-cortical fiber system is in the simplest case 
described by the following two-dimensional telegraph equation, 

(d \ 2 3 

— + v ek A ek \ <p ek (r, t) - -v 2 et y 2 <p ek (r, f) 

= v 2 ek A 2 ek N« k S e [h e (r,t)] (6) 
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where N" k is the total number of excitatory connections formed by 
long-range cortico-cortical axons on long on local population k, 
and assumes a single axonal conduction velocity v ek and an expo- 
nential fall off with distance (characteristic scale = UA ek ) of the 
strength of cortico-cortical connectivity. For simplicity, and given 
the fact that at least in anesthesia BS appears to have a degree of spa- 
tial uniformity, we chose to only study the spatially homogeneous 
case, i.e., v 2 = 0. 

Equations (l)-(4) and (6) represent a system of 8 coupled 
non-linear partial differential equations that typically defines the 
Liley model of electrocortical rhythmogenesis, which is capable of 
reproducing the main features of spontaneous human electroen- 
cephalogram (alpha resonance, "1//' activity). Table 1 summarizes 
all model parameters, their definitions, and approximate ranges. 

2.2. MODEL PARAMETERIZATION: GENERATION OF NORMATIVE 
PARAMETER SETS 

Because BS activity (at least that induced by anesthetic and seda- 
tive action) is assumed to ultimately arise out of a background of 
normal electroencephalographic activity it is important to define 
parametrically normative states in order to study how they may 
be perturbed during health and disease. We therefore chose to 



Table 1 | List of spatially averaged parameters for different types k, I s 
{e, /} of neuronal target populations in the electrocortical model of 
Liley et al. (Liley et al., 1999, 2002, 2011; Bojak and Liley, 2005; Frascoli 
et al., 2011), with typical ranges that are assumed to be 
physiologically admissible. 





Definition 


Min 


Max 


Units 




Resting membrane potential 


-80 


-60 


mV 


*t 


Passive membrane decay 


5 


150 


ms 




time 










Excitatory reversal potential 


-20 


10 


mV 


A? 


Inhibitory reversal potential 


-90 


h[-5 


mV 




EPSP peak amplitude 


0.1 


2.0 


mV 


r* 


IPSP peak amplitude 


0.1 


2.0 


mV 


VYek 


EPSP rise time to peak 


1 


10 


ms 


VVik 


IPSP rise time to peak 


2 


100 


ms 


N ek 


Number of excitatory 


1000 


5000 






cortico-cortical synapses 










Number of excitatory 


2000 


5000 






intracortical synapses 










Number of inhibitory 


100 


1000 






intracortical synapses 








Vek 


Axonal conduction velocity 


0.1 


1 


mm ms -1 




Decay scale cortico-cortical 


10 


100 


mm 




connectivity 








cmax 
^k 


Maximum firing rate 


0.05 


0.5 


ms -1 


P-k 


Mean firing threshold 


-55 


-40 


mV 


Ok 


Firing threshold standard 


2 


7 


mV 




deviation 








Pik 


Extra-cortical synaptic input 


0 


10 


ms -1 




rate 









utilize previously defined parameter sets (Bojak and Liley, 2005) 
that have the following properties: (i) are confined to the physio- 
logically admissible parameter space (see Table 1), (ii) give rise to 
electroencephalographically and physiologically plausible activity 
(" II f decay at low frequencies plus a relatively sharp peak at alpha 
frequencies, 8-13 Hz; mean excitatory/inhibitory neuronal firing 
rates <20s _1 ) and (hi) that exhibit transient increases in total 
power and monotonic reductions in mean frequency with respect 
to modeled anesthetic action (see below). In general such sets can 
be found by randomly searching the high dimensional physio- 
logically admissible (and plausible) parameter space. For further 
details see Bojak and Liley (2005). 

2.3. MODELING ANESTHETIC ACTION 

The range of molecular and cellular targets identified to date as 
sites of anesthetic action is so varied that a unitary biological mech- 
anism for anesthetic effect seems unlikely. Nevertheless, at least 
functionally, at the level of cortex anesthetics seem to act princi- 
pally by enhancing the actions of inhibitory activity (Liley et al, 
2011). Indeed from the perspective of the mean field model we 
have described many of its parameters can be related in a fairly 
straightforward way to these alterations in inhibitory activity and 
other identified sites of anesthetic action in cortex (see Table 2). 
However a parametrically more flexible description of the PSP, 
than is presently incorporated, is required to meaningfully model 
the effects of anesthetics in which we can independently vary peak 
amplitude, rise (t pe ak) and decay times. For example isoflurane, 
a volatile halogenated anesthetic, has been shown to prolong the 
decay time of the unitary IPSP without altering its time to peak. 
Fortunately a simple modification of the equation describing the 
dynamics of the PSP enables independent adjustment of the peak 
amplitude, rise (t pea k) an d decay times. By defining l\ k to satisfy 



3t 



+ Ylk(sik) 



dt 



+ Ylk(£lk) 



] Ilk(r, 



t) 



= Yik(sik) exp[j/ tt (e ft )/K,°]r tt A tt (r, f)> 
Ylk(eik) = £lkY? k /(e e,k ~ 1), Ylk = Ylk(s)e £,k 



(7) 
(8) 



Table 2 | Relationship between major experimentally identified sites 
of cortical anesthetic action and parameters of the electrocortical 
model of Liley et al. (Liley et al., 1999, 2002; Bojak and Liley, 2005). 



Site of action 


Main anesthetic effect 


Parameters 


2PK channels and 


Increase in tonic inhibition 


Pik. h r k 


extrasynaptic GABAa 






nACh receptors 


Reduction in tonic excitation 


Pek, h[ 


Synaptic GABA A 


Increase of IPSPs 


Yik. r ik 


AMPA/kainate receptors and 


Reduction of EPSPs 


Yek, ?ek 


NMDA receptors* 






Myelinated axons 


Slowdown of conduction* 


Vek 


Na channels 


Alteration of neuronal firing 





Table adapted from Liley et al. (2011). EPSR excitatory PSP; IPSP inhibitory PSP 



'Parameters will depend on membrane potential in this case. 

' Effect demonstrated in periphery speculative in cortex (Swindale, 2003). 
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where defines the time to peak, we can control the decay of 
the unitary PSP by altering > 0. Increasing will monoton- 
ically increase the decay time of the tail of the unitary PSP (see 
lower left panel, Figure 6). Empirically it is found that increas- 
ing the aqueous concentration of a range of GABAergic anesthetic 
agents leads to a progressive increase in the decay time of the uni- 
tary inhibitory PSP (e.g., Banks and Pearce, 1999) and thus 
will be a monotonic function of anesthetic concentration c, i.e., 
ejjfc(c). Liley et al. (2011), based on a range of empirical evidence, 
have numerically estimated £;j;(c) for the volatile anesthetic agent 
isoflurane. However because en{c) is not currently known for 
other GABAergic anesthetic agents we will assume that in general 

Ej'J-OC C. 

It is worth noting that equation (7) reduces to equation (2) as 
£ft— > 0. Further details regarding this formulation can be found 

in Bojak and Liley (2005). 

2.4. DEFINING A THEORETICAL BASIS FOR BURSTING 

We call a bursting system fast-slow if it can be written in the 
following form 

x = f(x,y) (fast oscillatory system) (9) 
y = fig(x,y) (slow modulatory system) (10) 

where x e R m describes the m-dimensional system responsible for 
the fast oscillatory (spiking in single neuron models) dynamics 
and yeR" the n-dimensional slow system that modulates the fast 
oscillations (or spiking behavior). The parameter \i represents the 
ratio of the time scales between the slow and fast system. It is 
typically assumed that Because fi can be made as small as 

we like equations (9) and (10) represent a singularly perturbed 
system. 

We will assume that equations (l)-(4) and (6), which puta- 
tively describe the genesis of the "fast oscillatory" resting EEG, 
correspond to the m-dimensional fast system. To motivate the 
slow n-dimensional slow modulatory system we will make a plau- 
sible appeal to the biophysics of synaptic resource depletion and 
recovery during periods of sustained neuronal population activity. 
During periods of high firing neural activity a variety of factors 
come into play to diminish synaptic efficiency. The most important 
of these are receptor desensitization and synaptic vesicle deple- 
tion. Tsodyks and Markram (1997) developed a model to account 
for the biophysics of such activity dependent short term synaptic 
depression estimating that its onset is rapid, of the order of mil- 
liseconds, but that its recovery is quite slow, of the order of 800 ms. 
Given that such a time scale is approximately at least an order of 
magnitude greater than the characteristic time scales associated 
with resting EEG activity, this may represent a candidate slow 
EEG modulatory system. On this basis we choose to include this 
activity dependent short term synaptic depression using the fol- 
lowing two different formulations, referred respectively to as SSI 
and SS2, 

fik = m\Qi-hSi(h)l (ssi) (n) 
r tt = wtr" i(\ + expMh, - m - r ft ], (SS2) (12) 



where 1/yU.; is the characteristic time scale of the respective slow 
modulatory system. Based on Tsodyks and Markram (1997) we 
will fix l//z; to 1000 ms. The advantage of the first formulation is 
that the rates of synaptic recovery (/x;f?/) and depletion (/x/fc/) can 
be separately specified. The advantage of the second formulation 
is that remains bounded between rSl (the resting value) and 
zero decreasing monotonically with increasing mean soma mem- 
brane potential hi, with rS remaining as a free parameter. For 
low levels of the respective neuronal activity (hi) there is very little 
decrease in V ^ until a threshold £ ; is reached, with the parame- 
ter k i controlling the sensitivity of the change at this threshold to 
variations in neuronal activity. This formulation has previously 
been used by Tabak and Rinzel (2005) in their mean field model 
for spontaneous electrical bursting activity in embryonic chick 
spinal cord. 

2.5. COMPUTATIONAL METHODS 

All numerical integrations and one-dimensional dynamical con- 
tinuations were performed using the XPPAUT package (Ermen- 
trout, 2002). A 4-th order Runge-Kutta scheme with a time step of 
0. 1 ms was used to numerically integrate the differential equations. 
Because of the multiple time scales our differential system may 
suffer from stiffness and numerical solutions may not converge. 
In these cases we have used the recommended "stiff" integra- 
tor CVODE (Cohen and Hindmarsh, 1994) as implemented in 
XPPAUT. 

3. RESULTS 

All numerical simulations were performed using a single model 
parameter set having a physiologically plausible white noise fluctu- 
ation spectrum (see top left panel Figure 3) and a single stable fixed 
point. This parameter set was found using the methods described 
in section 2.2. The parameters used, all within the physiologically 
admissible domain, can be found in the Table 3. 

Figure 3 shows the effects of the activity dependent modulation 
of r,-, on simulated mean field EEG activity using SSI [equation 
(11)]. Bursts emerge periodically, with intervening near isoelec- 
tric intervals, apparently driven by slow variations in r,-,- (red line, 
top right panel). A spectrogram of a sufficiently long simulated 
time series reveals that the frequency of the model EEG activity 
decreases from low beta (= 15 Hz) to high alpha (= 12 Hz) over 
the period of the bursts. Such intra-burst "chirping" is a common 
feature of many of the topologically identified single neuronal 
bursters (Izhikevich, 2007). Of interest are the multiple harmonics 
of this dominant oscillatory activity. 

Figure 4 illustrates how we might dynamically account for the 
transition to, and cessation from, bursting and follows the now 
standard method of the dissection of neural bursting pioneered by 
Rinzel ( 1985). Here we have set m = 0 and consider how the "fast" 
EEG system responds. Figure 4 shows a one-dimensional bifurca- 
tion diagram of this "fast" EEG system with r„ as the bifurcation 
parameter. Thick black lines show the fixed points as a function 
of r„. For small values of r,-,- there is a single stable fixed point. 
As Fa is increased this fixed point loses stability by a super-critical 
Hopf bifurcation, thus signaling the onset of limit cycle activity. 
Periodic continuations of this low amplitude activity reveals that 
it, and a stable fixed point, co-exist with a higher amplitude limit 
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FIGURE 3 | Effects of the slow, activity dependent modulation of 
IPSPs on the inhibitory neuronal population, r„ using the SS1 
model of equation (11 ). The top left panel shows the fluctuation 
spectrum of the unmodulated "fast" EEG system (Bojak and Liley, 
2005). When r, t of this system is slowly modulated (/i,- = 0.001 ms~' , 
fl, = 0.1818mV, k,= 10mVs) bursting emerges (top right panel). Bursting 



is associated with the periodic modulation of P, (red line, top right panel) 
on a much longer time scale than that of intra-burst oscillatory activity. 
Time-frequency analysis (bottom right panel) reveals that intra-burst 
activity sweeps down through a range of frequencies from low beta 
(13-30 Hz) to high alpha (8-12 Hz). Parameter values used for the "fast" 
EEG system can be found in Table 3. 



cycle, thus suggesting that an activity dependent hysteresis drives 
the system between a fixed point and a high amplitude oscilla- 
tion, thus giving rise to the bursting activity observed. This can be 
better seen by superimposing on this diagram the trajectory of a 
single burst (thin solid black line). Here we can see that the burst 
terminates through a fold-limit cycle bifurcation. At this stage it 
is not clear what bifurcation accounts for the emergence of the 
burst. 

A well described feature of anesthetic action is the reduc- 
tion in cerebral blood flow and metabolism (Kaisti et al, 
2003). Therefore during anesthetic action it would be reason- 
able to assume that the recovery of pre-synaptic neurotransmit- 
ter levels will be impaired. In particular as the anesthetic level 
increases then the rate of synaptic recovery should decrease. 
Figure 5 shows the effects of systematic reductions in the synap- 
tic recovery rate for the SSI model. As Q{ (synaptic recovery) 
is decreased the burst duration decreases and the period of 
non-oscillatory isoelectric activity increases until the model EEG 
becomes fully isoelectric. This trend is also observed clinically 
during anesthesia. 

Short term synaptic depression would be expected to affect all 
synapses, though the depression would not be expected to be uni- 
form. So far we have assumed that the synaptic depression would 
principally affect inhibitory synapses between inhibitory neu- 
rons. Will such bursting survive when all types of "fast" synaptic 
activity is subject to the biological forces of short term synap- 
tic depression? Figure 6 reveals that bursting does occur when 
both excitatory and inhibitory synaptic activity undergoes activity 
dependent short term synaptic depression. Further, the bursting 
that emerges is strongly modulated by parameters of the "fast" 



Table 3 | Model parameter set used in simulations of Figures 3-8. 



Parameter 


Value 


Parameter 


Value 


h'i st (mV) 


-68.1355 




4994.4860 


hf st (mV) 


-77.2602 




2222.9060 


hf e (mV) 


-15.8527 




4582.0661 


h e J (mV) 


7.4228 




4198.1829 


hf g q (mV) 


-85.9896 




989.5281 


hf (mV) 


-84.5363 




531.9419 


r e (ms) 


138.3660 


Vgk (cm ms 1 ) 


0.1714 


Tj (ms) 


89.3207 


A ek (cm -1 ) 


0.2433 


r ee (mV) 


0.3127 


S™ x (ms- 1 ) 


0.2801 


r e , (mV) 


0.9426 


S max (ms _l) 


0.1228 


r, e (mV) 


0.4947 


li e (mV) 


-47.1364 


r, v (mV) 


1.4122 


Hi (mV) 


-45.3751 


y ee (ms"') 


0.4393 


tie (mV) 


2.6120 


Yei (ms -1 ) 


0.2350 


a; (mV) 


2.8294 


Yie (ms" 1 ) 


0.0791 


Pee (ms~ 1 ) 


3.6032 


Ya (ms -1 ) 


0.0782 


p ej (ms -1 ) 


0.3639 



Parameter set taken from Bojak and Liley (2005). Because parameters were 
obtained as part of a numerical search their full precision had been detailed, 
however their sensitivity to perturbation is much less than the precision reported 

(Bojak and Liley, 2005). 



EEG system (see Table 2) that have been identified as targets 
for the action of anesthetic agents. Prolonging the decay of the 
unitary IPSP and reducing subcortical input are both found to 
significantly modulate modeled EEG bursting. In particular it is 
found that reducing excitatory extra-cortical input (p ee ), which 
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r.. (mV) 
II v ' 

FIGURE 4 | One-dimensional continuation of the "fast" EEG system in 
h e suggests how we might account for the dynamics of bursting. Thick 
black lines (solid = stable, dashed = unstable) shows the fixed points (fp) as 
a function of r„. Solid black circle indicates fixed point for the default 
parameter values for "fast" EEG system (parameters as for Figure 3). Solid 
blue lines indicate the locus of maximum and minimum amplitudes 
respectively of stable limit cycle (Ic) activity, whereas solid red lines 
correspond similarly to locus of unstable limit activity. Transitions between 
stable and unstable Ic activity is predominantly through fold-limit cycle 
bifurcations, except for middle limit cycle branch where transition is through 
a torus bifurcation. Superimposed on this one-dimensional bifurcation 
diagram is the trajectory of a single burst. 



presumably dominantly arises from thalamus, leads to very long 
quiescent (isoelectric) periods. 

One of the limitations in using SSI is that the PSPs are, in prin- 
ciple, free to take any value whereas physiology would dictate that 
they should remain bounded. To explore the effects of this restric- 
tion we chose to define an alternative slow modulating system [SS2, 
equation (12)]. Figure 7 shows the effect of utilizing this system 
to provide a slow activity dependent modulation of excitatory and 
inhibitory synaptic efficacy of our "fast" EEG system. In the left 
panel of this figure we have plotted c e = r e //r^ as a function of 
time. The interesting thing to note is that in contrast to Figure 3 
excitatory synaptic efficacy decreases during the quiescent inter- 
burst period and increases during the burst. The left hand panel 
however shows that there is a significant phase difference between 
the normalized excitatory synaptic efficacy c e =T e ;/r°, and the 



6 = 0.1818 

e = o.07 

1 ill 

e = 0.04 

2s 



FIGURE 5 | The effect of parametrically varying the rate of recovery of 
synaptic efficacy in modulating the modeled burst duration and the 
quiescent (isoelectric) period. If $, is made small enough the bursting 
solutions will undergo a bifurcation to a non-oscillatory state. All parameters 
as per Figure 4 



normalized inhibitory synaptic efficacy c, = r,-//r° ; . This suggests 
that there is an important dynamical interplay between excitatory 
and inhibitory synaptic efficacy to regulate neuronal population 
excitability such that bursting occurs. This may explain why there is 
confusion in the empirical literature regarding the role alterations 
in synaptic efficiency have in the genesis of BS. 

An important difference between SS2 and SS 1 is that parameters 
hypothesized to be important targets of anesthetic action p ee and 
£,-; are able not only to parametrically regulate bursting but appear 
also able to switch bursting on (presumably through a bifurcation 
from a stable fixed point). Figure 8 illustrates this. If the inhibitory 
neuronal IPSP decay time is not long enough then a single fixed 
point dominates which has an associated white noise fluctuation 
spectrum. But as the IPSP decay time increases (beyond eg > 1.8 
for the parameter set chosen) then bursting emerges. However 
if p ee is decreased, as we would expect during anesthesia, then 
the isoelectric period is prolonged until at some critical value of 
pee bursting is extinguished to be replaced by an infinitely long 
quiescent/isoelectric period. 

4. DISCUSSION 

We have described here how a well-known model of the "fast" 
dynamics of the EEG can be modulated by a number of slow sys- 
tems to produce bursting activity that bears some resemblance to 
BS seen clinically. The slow systems that we used were all based 
on some form of activity dependent short term plasticity that has 
been empirically observed, and used successfully in other models 
of macroscopic level bursting (Tabak and Rinzel, 2005). While we 
were able to clearly show the existence of bursting, because we 
did not include any additive or multiplicative noise sources, we 
were unable to account for the quasi-periodicity of BS. Thus all 
our busting arises from purely deterministic processes, presum- 
ably involving a range of well described bifurcations (Izhikevich, 
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FIGURE 6 | The effects of short term synaptic depression of both 
excitatory and inhibitory cortical synapses in the genesis of burst 
suppression and its modulation by variations in extra-cortical input (p ee ) 



and the IPSP decay time [s n , s ie ) (see section 2.3 for further details). 

Parameters: e„ = 0.1818 mV, e, = 0.07 mV, k e = 14mVs. All other parameters 
as for Figure 3. 




- 0.85 


0.85 


0.82 




■ 0.79 


0.8 






■ 0.76 


0.75 


0.73 




- 0.7 


07 



time 




0.95 



FIGURE 7 | Bursting produced by the activity dependent slow 
modulation of excitatory and inhibitory synaptic efficacy 
according to SS2 [equation (12)]. Left panel shows h, (solid black 
line) and normalized excitatory synaptic efficacy c B = r s ,/r°, (solid red 



line) as a function of time. The left hand panel shows the phase 
relationship between normalized excitatory and inhibitory (c, = r,/r,°) 
synaptic efficacies. Parameters: k, = 0.2 mV~\ k, = 0.1 mV~\ e, = 1.8, 
Si, = 1.5. 
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FIGURE 8 | Parametric variability of bursting produced by the activity dependent slow modulation of excitatory and inhibitory synaptic efficacy 
according to SS2 [equation (12)]. For s„ < 1.8 simulated EEG was isoelectric as was the case for labeled combinations of s, and p„. Parameters as for Figure 7 



2007). However because our system clearly exhibited bistability 
(see Figure 4) it is almost certain that our system will be able to 
exhibit some form of burst excitability in response to stochastic 
forcing. Such burst excitability has been described experimentally. 
For example during BS induced by various halogenated anes- 
thetic agents, bursts can be readily evoked by auditory, visual, 
or somatosensory stimuli (Hartikainen et al., 1995). Neverthe- 
less while burst onset and duration may be random variables it 
would seem that the bursts themselves should reveal a high degree 
of determinism (weak non-linearity) when compared to EEG in 
which bursting or epileptiform activity is not evident. 

Because the parameter space of the underlying "fast" EEG 
model is potentially extremely large it is not possible to system- 
atically explore its dynamical repertoire and it may be possible 
that this system, not augmented with one of the slow systems 
described, is able to burst. Nevertheless, on the basis of our 
results, and what is known regarding the dynamical mechanisms 



of bursting, it would seem likely that multiple pathways to BS 
exist through a variety of activity dependent slow modulatory 
systems. 

Further we might hypothesize that such slow modulatory sys- 
tems might span a number of functional scales in the brain. 
Figure 9 diagrammatically illustrates some possible candidate sys- 
tems. An obvious activity driven slow modulatory system would 
be that associated with thalamus and the corresponding thalamo- 
cortical feed-back. Mean field models of the EEG that incorporate 
thalamo-cortical feed-back have been developed (Rennie et al., 
2002) and it will be interesting to see if they are structurally config- 
ured to support BS. In addition to synaptic fatigue another obvious 
cortical system that might be marshaled to provide slow activity 
dependent modulation is the cortico-cortical conduction system. 
Although at this stage there is currently little evidence to suggest 
anesthetics slow conduction velocities, it is widely documented 
that axonal conduction velocities are significantly decreased in 
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FIGURE 9 | Diagram illustrating a number of hypothesized, 
physiologically plausible, multiple scale slow modulatory systems 
that could be important in the genesis of electroencephalographic 
bursting in the context mean field model of Liley et al. (A) 

Thalamo-cortical feed-back (B) changes in the conduction properties of the 
long-range cortico-cortical fiber system (C) slow changes in the efficacy of 
synaptic function due to activity dependent resource depletion and 
restitution. See accompanying text for symbol definitions. 



hypothermia. Decreases of up to 5% per °C for conduction veloc- 
ity have been reported (Waxman, 1980). As BS has been observed 
to occur in hypothermia (Schwartz et al., 1989; Akrawi et al., 1996) 
we can conclude that a slow system emerging in the long-range 



coupling via a slowing of axonal conduction velocity is a possible 
route to BS. 

Our attempts to account for the dynamical pathogenesis of 
burst suppression differ from other approaches, most notably 
Ching et al. (20 12 ) . Ching et al. described the scalp EEC in terms of 
the activity of a small scale, biophysically detailed, computational 
model of interacting populations of cortical and thalamic neurons. 
Burst suppression was modeled as arising from metabolically 
induced alterations in an ATP-gated slow neuronal membrane 
potassium current (Ik atp ) on the basis that the reduction in cere- 
bral metabolic rate (CMR02) induced by anesthetic agents and 
hypoxia was associated with the depletion of ATP, and hence mem- 
brane hyperpolarization. While on this basis they claim to have 
accounted for a number of defining features of BS that included 
(i) the spatial synchrony of burst onset (ii) the parametric variabil- 
ity of burst duration/isoelectricity and (iii) the characteristically 
long time scales associated with bursting/isoelectricity compared 
to resting EEG, some caution needs to be exercised. 

Firstly their model of resting/spontaneous EEG is constructed 
on the basis of the activity of no more than 20 model neurons. 
Because EEG is a distributed large scale phenomena such a model 
is unlikely to meaningfully account for resting/spontaneous activ- 
ity particularly given the absence of any long-range excitatory 
cortico-cortical coupling. This has important implications for the 
propagation of burst activity particularly given that the onset of 
bursts, when examined at fine temporal scales, is probably not 
truly spatially homogeneous. 

Secondly while the relationship between CMR02 and ATP pro- 
duction cannot be reasonably disputed, not all anesthetic agents 
that produce reductions in CMR02 produce BS. For example the 
noble gas xenon has been reported to reduce CMR02 by up to 
33% in human participants (Rex et al, 2006) yet is not associated 
with any discernable BS. 

Thirdly the approach they have taken to producing BS essen- 
tially depends on the slow modulation of a faster system, the 
approach we have adopted here. The modeled time scales of Ik A t P 
variability are very long, of the order of tens of seconds. 

For clarity and tractability the current investigations have 
focused on the spatially homogeneous case for the model of Liley 
et al. [i.e., y 2 = 0 in equation (6)]. Clearly the emergence of BS in 
the spatially extended case will need to be investigated through the 
appropriate numerical solution of the defining partial differen- 
tial equations. Because the cortical phase synchrony (Hartikainen 
et al, 1995) of burst suppression has not, as far as we are aware, 
been explicitly investigated it will be crucial to empirically deter- 
mine the spatiotemporal emergence of bursts in order to assess 
the importance of excitatory cortico-cortical connectivity in the 
emergence and modulation of BS as implied by the model and as 
we have suggested. This will require recording high density EEG 
during anesthesia in which BS is present. 
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